{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://habrastorage.org/webt/3x/l-/yv/3xl-yvgopmrwrjaz7hvkouaplye.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Автор: Лопухова Арина (Slack: @erynn)_\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Содержание <a name=\"cont\"></a>\n",
    "1. [Описание задачи и данных](#intro)\n",
    "2. [Создание признаков I](#feat_1)\n",
    "3. [Первичный и визуальный анализ. Найденные зависимости](#analysis)\n",
    "4. [Создание признаков II](#feat_2)\n",
    "5. [Выбор моделей и метрики](#model)\n",
    "6. [Предобработка данных, обучение, оценка результатов на LB](#learning)\n",
    "    \n",
    "    6.1 [KNeighborsRegressor](#knr)\n",
    "    \n",
    "    6.2 [KNeighborsRegressor на расширенной выборке](#knr)\n",
    "    \n",
    "    6.3 [RandomForestRegressor](#rfr)\n",
    "    \n",
    "    6.4 [XGBRegressor](#xgb)\n",
    "\n",
    "7. [Почти выводы](#alm_conc)\n",
    "8. [Выводы](#conc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Описание задачи и данных <a name=\"intro\"></a>\n",
    "\n",
    "[К началу](#cont)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Подробное описание задачи и данных можно посмотреть [на странице соревнования.](https://datasouls.com/c/mts-geohack/description) Задача в том, чтобы по среднему числу звонков в экстренные службы с участков примерно 500x500 метров в западной части Москвы и Московской области спрогнозировать число звонков с таких участков в восточной части Москвы и Московской области по каждому дню недели. Данные 2017 года.  \n",
    "\n",
    "Качество оценивается через [коэффициент ранговой корреляции Кендалла](http://www.machinelearning.ru/wiki/index.php?title=%D0%9A%D0%BE%D1%8D%D1%84%D1%84%D0%B8%D1%86%D0%B8%D0%B5%D0%BD%D1%82_%D0%BA%D0%BE%D1%80%D1%80%D0%B5%D0%BB%D1%8F%D1%86%D0%B8%D0%B8_%D0%9A%D0%B5%D0%BD%D0%B4%D0%B5%D0%BB%D0%BB%D0%B0) независимо по каждому дню недели. \n",
    "\n",
    "Задача интересна тем, что потенциально может выявить проблемные районы, которым нужны изменения в инфраструктуре или больше внимания экстренных служб. Она также полезна как анализ нагрузки на сотового оператора (все-таки данные предоставляет МТС). Плюс задача хороша в учебно-исследовательском плане, потому что предполагает работу с геоданными, к тому же из внешних источников вроде [OpenStreetMap](https://www.openstreetmap.org/). Готовых признаков вообще по сути нет, и все нужно делать вручную. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Данные (train + test в одном файле и все остальные дополнительные данные) можно скачать [отсюда](https://drive.google.com/file/d/10F166_vlWkP2-z526nTaiWG4efMyIf8a/view?usp=sharing).\n",
    "\n",
    "Значения переменных:\n",
    "* __lat_bl, lon_bl__ — координаты нижнего левого угла участка\n",
    "* __lat_tr, lon_tr__ — координаты верхнего правого угла участка\n",
    "* __lat_c, lon_c__ — координаты центра участка\n",
    "* __calls_daily__ — среднее число вызовов по всем дням\n",
    "* __calls_workday__ — среднее число вызовов по рабочим дням\n",
    "* __calls_weekend__ — среднее число вызовов по выходным дням\n",
    "* __calls_wd{D}__ — среднее число вызовов по дню недели D (0 — понедельник, 6 — воскресенье)\n",
    "* __is_test__ — индикатор того, что объект тестовый (у тестовых объектов значения числа звонков не показаны)\n",
    "* __is_target__ — индикатор целевого участка (указан только для обучающей выборки)\n",
    "\n",
    "Про целевой участок: качество предсказания оценивается не по всей тестовой части, а только по тем участкам, которые входят в целевые. Участки из тестовой выборки отбираются в целевые по тому же принципу, что и участки обучающей (по какому — участникам неизвестно).\n",
    "\n",
    "Целевой признак — __calls_wd{D}__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Примечание к работе:__ в [репозитории](https://github.com/datasouls/mts-geohack/blob/master/Geohack112_StarterKit.ipynb) соревнования есть базовая тетрадка Starter Kit, части которой здесь будут использоваться. Оттуда взяты идеи про расстояние до Кремля и близость к каким-то объектам/число объектов в радиусе (часть про конструирование признаков). Взятый из Kit код будет помечен как #SK."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Установка osmread для работы с OpenStreetMap и импорт библиотек__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install osmread"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import pickle\n",
    "from sklearn.neighbors import NearestNeighbors\n",
    "from sklearn.model_selection import train_test_split\n",
    "from scipy.stats import kendalltau\n",
    "import osmread\n",
    "\n",
    "from tqdm import tqdm_notebook\n",
    "import seaborn as sns\n",
    "import matplotlib as mpl\n",
    "from matplotlib import pyplot as plt\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#размер графиков\n",
    "mpl.rcParams['figure.figsize'] = (8, 7)\n",
    "#стиль и размер шрифта\n",
    "sns.set(style='whitegrid', font_scale=1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PATH_TO_DATA = './data/'  #путь к данным "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "distance из cython_dist — реализованная на cython функция для подсчета расстояния в км между точками в географических координатах"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cython_dist import distance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Считаем файл с участками__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zones_df = pd.read_csv('data/zones.csv', index_col='zone_id')\n",
    "zones_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Считаем из tagged_nodes.pickle объекты из OSM__, которые попали в интересующий нас регион и имеют хотя бы один тег. Подробно о том, как шла работа с картой и создавался tagged_nodes.pickle, можно посмотреть в [Starter Kit](https://github.com/datasouls/mts-geohack/blob/master/Geohack112_StarterKit.ipynb), раздел \"Работа с данными OpenStreetMap\".  \n",
    "\n",
    "Если вкратце про объекты карты OSM, то их три типа (описание из того же Starter Kit): \n",
    "* __Node__  — точки на карте\n",
    "* __Way__  — дороги, площади, задающиеся набором точек\n",
    "* __Relation__  — связи между элементами\n",
    "\n",
    "В tagged_nodes.pickle только объекты типа Node, которые попали в большой прямоугольник от Наро-Фоминска на юго-западе от Москвы до Красноармейска на северо-востоке от Москвы."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.path.join(PATH_TO_DATA, 'tagged_nodes.pickle'), 'rb') as f:\n",
    "    tagged_nodes = pickle.load(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Пример объекта__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tagged_nodes[17]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Создание признаков I <a name=\"feat_1\"></a>\n",
    "\n",
    "[К началу](#cont)\n",
    "\n",
    "С этого все начнется, потому что использовать готовые данные в виде чистых координат не имеет смысла. Уже после накидывания признаков будет проведен их анализ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Расстояние до Кремля__\n",
    "\n",
    "Мотивация: районы ближе к центру воспринимаются как более благополучные по сравнению с окраинами, что должно влиять на количество звонков в полицию.\n",
    "\n",
    "__Признаки вида: 1) число объектов какого-то типа в определенном радиусе от центра участка; 2) расстояние до k ближайших объектов какого-то типа; 3) расстояние от центра участка до ближайшего объекта определенного типа__\n",
    "\n",
    "Под объектами имеются в виду объекты OpenStreetMap.\n",
    "\n",
    "Мотивация:\n",
    "* __Магазины__: 1) плотная застройка магазинами позволяет отделить более старые обжитые районы от новых микрорайонов с новостройками ближе к окраинам; 2) магазины — открытые многолюдные места, которые обеспечивают поток людей на ближайших улицах, что должно уменьшать вероятность преступлений на улице \n",
    "* __Рестораны, кафе + кино, театр + каток, площадка, парк__: аналогичный поток людей и признак благоустроенности\n",
    "* __Общественный транспорт__: вот тут сложнее, потому что с одной стороны это близость людей, с другой — особенно плотное скопление людей, давка и потенциальная опасность здоровью плюс возможность ДТП, поэтому пока неясно, в какую сторону этот признак может влиять на целевой"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Признаки могут считаться довольно долго из-за пользовательской метрики расстояния, которая передается в NearestNeighbours, поэтому они уже заранее посчитаны и лежат в dist_features.pkl. Про расшифровку тегов объектов OSM можно почитать\n",
    "[здесь.](https://wiki.openstreetmap.org/wiki/Map_Features)\n",
    "\n",
    "Примечание про метрику расстояния: в исходном SK расстояние между точками в градусах считалось как евклидово, в моей версии оно считается именно как расстояние между точками на эллипсоиде ради интерпретируемости и возможного улучшения качества."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_features = pd.DataFrame(index=zones_df.index)\n",
    "\n",
    "#признак вида \"расстояние до Кремля\"\n",
    "kremlin_lat, kremlin_lon = 55.753722, 37.620657\n",
    "dist_features['kremlin_dist'] = zones_df.apply(lambda row: \\\n",
    "                                             distance(np.array([row.lat_c, row.lon_c]), np.array([kremlin_lat, kremlin_lon])), \\\n",
    "                                             axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SK\n",
    "#набор фильтров точек, по которым будет считаться статистика\n",
    "\n",
    "POINT_FEATURE_FILTERS = [\n",
    "    ('tagged', lambda node: len(node.tags) > 0),\n",
    "    ('railway', lambda node: node.tags.get('railway') == 'station'),\n",
    "    ('public_transport', lambda node: 'public_transport' in node.tags),\n",
    "    ('entertainment', lambda node: node.tags.get('amenity') in ['cinema', 'theatre']),\n",
    "    ('food', lambda node: node.tags.get('amenity') in ['restaurant', 'cafe', 'fast_food']),\n",
    "    ('leisure', lambda node: node.tags.get('leisure') not in [None, 'amusement_arcade', 'adult_gaming_centre']),\n",
    "    ('shop', lambda node: node.tags.get('shop') not in [None, 'lottery'])\n",
    "]\n",
    "\n",
    "# центры квадратов в виде матрицы\n",
    "X_zone_centers = zones_df[['lat_c', 'lon_c']].as_matrix()\n",
    "\n",
    "for prefix, point_filter in tqdm_notebook(POINT_FEATURE_FILTERS):\n",
    "\n",
    "    # берем подмножество точек в соответствии с фильтром\n",
    "    coords = np.array([\n",
    "        [node.lat, node.lon]\n",
    "        for node in tagged_nodes\n",
    "        if point_filter(node)\n",
    "    ])\n",
    "\n",
    "    # строим структуру данных для быстрого поиска точек\n",
    "    neighbors = NearestNeighbors(metric=distance).fit(coords)\n",
    "    \n",
    "    # признак вида \"количество точек в радиусе R от центра квадрата\"\n",
    "    for radius in tqdm_notebook([0.005, 0.01, 0.05, 0.1, 0.15, 0.2, 0.25]):\n",
    "        dists, inds = neighbors.radius_neighbors(X=X_zone_centers, radius=radius)\n",
    "        dist_features['{}_points_in_{}'.format(prefix, radius)] = np.array([len(x) for x in inds])\n",
    "\n",
    "    # признак вида \"расстояние до ближайших K точек\"\n",
    "    for n_neighbors in tqdm_notebook([3, 5, 10]):\n",
    "        dists, inds = neighbors.kneighbors(X=X_zone_centers, n_neighbors=n_neighbors)\n",
    "        dist_features['{}_mean_dist_k_{}'.format(prefix, n_neighbors)] = dists.mean(axis=1)\n",
    "        dist_features['{}_max_dist_k_{}'.format(prefix, n_neighbors)] = dists.max(axis=1)\n",
    "        dist_features['{}_std_dist_k_{}'.format(prefix, n_neighbors)] = dists.std(axis=1)\n",
    "\n",
    "    # признак вида \"расстояние до ближайшей точки\"\n",
    "    dist_features['{}_min'.format(prefix)] = dists.min(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Считаем признаки из файла__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.path.join(PATH_TO_DATA, 'features/dist_features.pkl'), 'rb') as f:\n",
    "    dist_features = pickle.load(f)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_features.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Переходы через проезжую часть/железную дорогу__:\n",
    "\n",
    "* число пешеходных переходов\n",
    "* число и доля нерегулируемых пешеходных переходов\n",
    "\n",
    "Мотивация: переход через дорогу, особенно в нерегулируемых местах, связан с риском аварии, что в свою очередь связано со звонками в скорую. Проверить, зайдет ли доля объектов в качестве признака, в голову пришло после тренировки по машинному обучению в Яндексе, и там это просто сработало в одном соревновании."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def in_zone(zone_coord, obj_coord):\n",
    "    '''\n",
    "    in: iterable zone_coord (lat_bl, lon_bl, lat_tr, lon_tr) \n",
    "        iterable obj_coord (lat, lon)\n",
    "    out: True, если объект внутри зоны, иначе False\n",
    "    '''\n",
    "    lat_bl, lon_bl, lat_tr, lon_tr = zone_coord\n",
    "    lat, lon = obj_coord\n",
    "    \n",
    "    return (lat_bl <= lat <= lat_tr) and (lon_bl <= lon <= lon_tr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def obj_count(zone_coord, nodes):\n",
    "    '''\n",
    "    in: iterable zone_coord (lat_bl, lon_bl, lat_tr, lon_tr)\n",
    "        list  nodes (список координат объектов одного интересующего типа)\n",
    "    out: число объектов указанного типа внутри участка\n",
    "    '''\n",
    "    count = 0\n",
    "    for node in nodes:\n",
    "        if in_zone(zone_coord, node):\n",
    "            count += 1\n",
    "    return count        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def select_type(types, tagged_nodes=tagged_nodes):\n",
    "    '''\n",
    "    in: str types (строка вида 'key1=value1,key2=value2, ...')  с интересующими парами ключ-значение\n",
    "        либо только ключом (тогда проверяется, что ключ есть среди ключей объекта OSM)\n",
    "        list tagged_nodes (объекты OSM)\n",
    "    out: list (список списков, где во внутренних списках — координаты объектов OSM, соответствующих парам из types)\n",
    "    '''    \n",
    "    outer, inner = [], []\n",
    "    for pair in types.split(sep=','):\n",
    "        try:\n",
    "            key, value = pair.split(sep='=')\n",
    "        except Exception:\n",
    "            key = pair\n",
    "            for node in tagged_nodes:\n",
    "                if key in node.tags:\n",
    "                    inner.append((node.lat, node.lon)) \n",
    "        else:            \n",
    "            for node in tagged_nodes:\n",
    "                if node.tags.get(key) == value:\n",
    "                    inner.append((node.lat, node.lon))\n",
    "\n",
    "        outer.append(inner.copy())\n",
    "        inner.clear()\n",
    "        \n",
    "    return outer    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crossing_nodes, crossing_uncontr_nodes = select_type('crossing,crossing=uncontrolled') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "n_crossing_nodes = zones_df.apply(lambda row: obj_count((row.lat_bl, row.lon_bl, row.lat_tr, row.lon_tr),\n",
    "                                                        crossing_nodes), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "n_crossing_uncontr_nodes = zones_df.apply(lambda row: obj_count((row.lat_bl, row.lon_bl, row.lat_tr, row.lon_tr),\n",
    "                                                      crossing_uncontr_nodes), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prt_crossing_uncontr_nodes = (n_crossing_uncontr_nodes / n_crossing_nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crossing_features = pd.DataFrame(data={'n_crossing': n_crossing_nodes,\n",
    "                                       'n_crossing_uncontr': n_crossing_uncontr_nodes,\n",
    "                                       'prt_crossing': prt_crossing_uncontr_nodes})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Можно считать из файла__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.path.join(PATH_TO_DATA, 'features/crossing_features.pkl'), 'rb') as f:\n",
    "    crossing_features = pickle.load(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Признаки, связанные с жилыми домами:__\n",
    "\n",
    "* число жилых домов в участке\n",
    "* число и доля жилых домов, введенных в эксплуатацию раньше 1970\n",
    "* число и доля жилых домов по индивидуальному проекту\n",
    "* число и доля жилых домов, к которым не проведен газ\n",
    "* число и доля жилых домов с автоматической пожарной системой \n",
    "\n",
    "Мотивация: \n",
    "\n",
    "* в старых домах обычно живет много пожилых людей, которым чаще может потребоваться скорая помощь (вообще тут можно засомневаться и вспомнить про реновацию в Москве, но в 2017 году массового переселения в новые дома еще не было)\n",
    "* индивидуальные проекты как характеристика района и, предположительно, благоустроенности \n",
    "* доступ к газу связан с риском ЧС\n",
    "* автоматические системы тушения, предположительно, могут влиять на частоту вызовов в пожарные службы\n",
    "\n",
    "Признаки получаются на основе данных из house_moscow_info.csv и house_nonmoscow_info.csv, которые были вытащены [отсюда](https://www.reformagkh.ru/opendata) (реестр домов по городу Москве и Московской области)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Считываем данные__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_moscow_info = pd.read_csv(os.path.join(PATH_TO_DATA, 'house_moscow_info.csv'),\n",
    "                         sep=';', \n",
    "                         usecols=['id', 'address', 'exploitation_start_year', \n",
    "                                  'project_type', 'gas_type', 'firefighting_type'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_moscow_info.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_nonmoscow_info = pd.read_csv(os.path.join(PATH_TO_DATA, 'house_nonmoscow_info.csv'),\n",
    "                         sep=';',\n",
    "                         usecols=['id', 'formalname_city', 'address', 'exploitation_start_year', \n",
    "                                  'project_type', 'gas_type', 'firefighting_type'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Из Московской области возьмем только топ-50 крупнейших городов: если посмотреть на целевые зоны в обучающей выборке и на карту, где эти зоны отмечены, то можно предположить, что целевыми считаются крупнейшие поселения. \n",
    "\n",
    "__Считаем список крупнейших городов МО__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.path.join(PATH_TO_DATA, 'top_c.txt')) as f:\n",
    "    top_c = f.readlines()\n",
    "    \n",
    "    \n",
    "top_c = [s.strip('\\n').strip('\\ufeff') for s in top_c] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_nonmoscow_info = house_nonmoscow_info[house_nonmoscow_info.formalname_city.isin(top_c)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_nonmoscow_info.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Приведем записи в текстовых колонках к одному формату__\n",
    "\n",
    "Например, значения вроде \"индивидуальный\", \"инд.\" и подобные в колонке project_type все переделаем в \"инд\", чтобы потом правильно отбирать дома по условию \"индивидуальный проект/не индивидуальный проект\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ind_proj_names_moscow = [s for s in list(house_moscow_info.project_type.value_counts().index) if 'инд' in s]\n",
    "ind_proj_names_moscow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_moscow_info.loc[house_moscow_info.project_type.isin(ind_proj_names_moscow), 'project_type'] = 'инд'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ind_proj_names_nonmoscow = [s for s in list(house_nonmoscow_info.project_type.value_counts().index) if 'инд' in s]\n",
    "ind_proj_names_nonmoscow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_nonmoscow_info.loc[house_nonmoscow_info.project_type.isin(ind_proj_names_nonmoscow), 'project_type'] = 'инд'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "У типов системы газоснабжения и противопожарной системы все в одном формате."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_moscow_info['gas_type'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_moscow_info['firefighting_type'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_nonmoscow_info['gas_type'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_nonmoscow_info['firefighting_type'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Сразу несколько проблем__\n",
    "\n",
    "У OSM проблемы с жилыми домами (дома не отмечают как жилые либо вообще нет тега \"здание\"), поэтому использовалась библиотека __geocoder__, которая связывалась с API Яндекс карт и получала координаты адресов из house_moscow_info и house_nonmoscow_info, чтобы на их основе создать новые признаки (координат в чистом виде среди признаков нет). \n",
    "\n",
    "Проблема в том, что выкладывать координаты в открытый доступ в виде файла нарушает правила пользования сервисом, поэтому дальше будет закомментированный код-пример того, как на основе координат создавались признаки. При этом самих координат в файлах проекта нет и вручную выполнить код в этих ячейках не получится. Готовый датафрейм с этими признаками можно считать из файла чуть ниже. \n",
    "\n",
    "Это все довольно чудовищно выглядит, но пока это лучшее решение проблемы, до которого удалось дойти.\n",
    "\n",
    "В house_moscow_lat_lon и house_nonmoscow_lat_lon лежали координаты адресов домов из house_moscow_info и house_nonmoscow_info."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# house_moscow_info['lat_lon'] = house_moscow_lat_lon\n",
    "# house_nonmoscow_info['lat_lon'] = house_nonmoscow_lat_lon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# full_house_info = pd.concat([house_moscow_info, house_nonmoscow_info.drop('formalname_city', axis=1)], \n",
    "#                             ignore_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Генерируем признаки__\n",
    "\n",
    "Попробуем только для целевых и тестовых квадратов. Получившиеся признаки можно сразу считать из house_features.pkl."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Среди координат в full_house_info попадаются -1: это когда Яндекс-карты не распознавали адрес. Такие адреса выкинем."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mask = full_house_info['lat_lon'].apply(lambda lst: type(lst) == int)\n",
    "# corrupted_coord = full_house_info['lat_lon'][mask]\n",
    "# full_house_info.drop(corrupted_coord.index, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%time\n",
    "# n_resident_house = zones_df[(zones_df.is_target == 1) | (zones_df.is_test == 1)] \\\n",
    "#                   .apply(lambda row: obj_count((row.lat_bl, row.lon_bl, row.lat_tr, row.lon_tr),\n",
    "#                                                 full_house_info['lat_lon']), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%time\n",
    "# n_old = zones_df[(zones_df.is_target == 1) | (zones_df.is_test == 1)] \\\n",
    "#         .apply(lambda row: obj_count((row.lat_bl, row.lon_bl, row.lat_tr, row.lon_tr),\n",
    "#                                       full_house_info[full_house_info.exploitation_start_year <= 1970]['lat_lon']), axis=1)\n",
    "# prt_old = n_old / n_resident_house"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%time\n",
    "# n_ind_proj = zones_df[(zones_df.is_target == 1) | (zones_df.is_test == 1)] \\\n",
    "#             .apply(lambda row: obj_count((row.lat_bl, row.lon_bl, row.lat_tr, row.lon_tr),\n",
    "#                                          full_house_info[full_house_info.project_type == 'инд']['lat_lon']), axis=1)\n",
    "# prt_ind_proj = n_ind_proj / n_resident_house    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%time\n",
    "# n_without_gas = zones_df[(zones_df.is_target == 1) | (zones_df.is_test == 1)] \\\n",
    "#                 .apply(lambda row: obj_count((row.lat_bl, row.lon_bl, row.lat_tr, row.lon_tr),\n",
    "#                                              full_house_info[full_house_info.gas_type == 'Отсутствует']['lat_lon']), axis=1)\n",
    "# prt_without_gas = n_without_gas / n_resident_house    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%time\n",
    "# n_auto_fire = zones_df[(zones_df.is_target == 1) | (zones_df.is_test == 1)] \\\n",
    "#              .apply(lambda row: obj_count((row.lat_bl, row.lon_bl, row.lat_tr, row.lon_tr),\n",
    "#                                           full_house_info[full_house_info.firefighting_type == 'Автоматическая']['lat_lon']), axis=1)\n",
    "# prt_auto_fire = n_auto_fire / n_resident_house    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# full_house_features = pd.DataFrame({'n_resident': n_resident_house, 'n_old': n_old, 'n_ind_proj': n_ind_proj,\n",
    "#                                     'n_without_gas': n_without_gas, 'n_auto_fire': n_auto_fire,\n",
    "#                                     'prt_old': prt_old, 'prt_ind_proj': prt_ind_proj,\n",
    "#                                     'prt_without_gas': prt_without_gas, 'prt_auto_fire': prt_auto_fire})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Тоже можно считать из файла__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.path.join(PATH_TO_DATA, 'features/house_features.pkl'), 'rb') as f:\n",
    "    full_house_features = pickle.load(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Число базовых станций МТС в участке__\n",
    "\n",
    "Координаты БС взяты [отсюда](http://bsmaps.ru/viewtopic.php?f=5&t=358) (оставлены только координаты нужного нам региона).\n",
    "\n",
    "Мотивация: предположение о том, что больше вышек там, где плотность звонков выше либо область значима."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.path.join(PATH_TO_DATA, 'stat_coords.pkl'), 'rb') as f:\n",
    "    stat_coords = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "n_base_stations = zones_df[(zones_df.is_target == 1) | (zones_df.is_test == 1)] \\\n",
    "                  .apply(lambda row: obj_count((row.lat_bl, row.lon_bl, row.lat_tr, row.lon_tr),\n",
    "                                                stat_coords), axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__И это можно считать из файла__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.path.join(PATH_TO_DATA, 'features/base_stat.pkl'), 'rb') as f:\n",
    "    n_base_stations = pickle.load(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Признак: есть ли в участке объект определенного типа__\n",
    "\n",
    "Мотивация: просто попробовать, потому что объекты значимые."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def obj_present(zone_coord, nodes):\n",
    "    '''\n",
    "    in: iterable zone_coord (lat_bl, lon_bl, lat_tr, lon_tr)\n",
    "        list  nodes (список координат объектов одного интересующего типа)\n",
    "    out: 1/0 (объект есть в участке/нет)\n",
    "    '''\n",
    "    for node in nodes:\n",
    "        if in_zone(zone_coord, node):\n",
    "            return 1\n",
    "    return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clinic_nodes, police_nodes, school_nodes = select_type('amenity=clinic,amenity=police,amenity=school') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "has_clinic = zones_df[(zones_df.is_target == 1) | (zones_df.is_test == 1)] \\\n",
    "            .apply(lambda row: obj_present((row.lat_bl, row.lon_bl, row.lat_tr, row.lon_tr),\n",
    "                                            clinic_nodes), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "has_police = zones_df[(zones_df.is_target == 1) | (zones_df.is_test == 1)] \\\n",
    "            .apply(lambda row: obj_present((row.lat_bl, row.lon_bl, row.lat_tr, row.lon_tr),\n",
    "                                            police_nodes), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "has_school = zones_df[(zones_df.is_target == 1) | (zones_df.is_test == 1)] \\\n",
    "            .apply(lambda row: obj_present((row.lat_bl, row.lon_bl, row.lat_tr, row.lon_tr),\n",
    "                                            school_nodes), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "has_obj_features = pd.DataFrame(data={'has_clinic': has_clinic,\n",
    "                                      'has_police': has_police,\n",
    "                                      'has_school': has_school})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Первичный + визуальный анализ данных\n",
    "\n",
    "## Найденные зависимости <a name=\"analysis\"></a>\n",
    "\n",
    "[К началу](#cont)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Создадим сводный DataFrame признаков и выделим train, test и звонки в train__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_full = zones_df[(zones_df.is_test == 1) | (zones_df.is_target == 1)].index\n",
    "full_features = pd.concat([dist_features.loc[idx_full],\n",
    "                           full_house_features,\n",
    "                           crossing_features.loc[idx_full],\n",
    "                           n_base_stations.rename('n_base_stations'),\n",
    "                           has_obj_features],\n",
    "                           axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_train, idx_test = zones_df[zones_df.is_target == 1].index, zones_df[zones_df.is_test == 1].index\n",
    "train_features, calls_train = full_features.loc[idx_train], zones_df.loc[idx_train, :].iloc[:, 8:] \n",
    "test_features = full_features.loc[idx_test]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_features.shape, test_features.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Посмотрим на значения NaN и незаполненные значения__\n",
    "\n",
    "Все признаки создавались вручную, так что можно сразу сказать, где есть проблемы. Незаполненных значений нет. NaN есть в тех колонках, которые отвечают за долю каких-то объектов. Появляются, очевидно, когда считается доля от нуля. Например:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_features['prt_old'].value_counts(dropna=False).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Диапазоны значений оценим по каждому типу признаков отдельно__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__По признакам расстояния__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_features.loc[idx_full].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_features.loc[idx_full].describe().iloc[:, 1:17]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Максимальное значение tagged_points напрягает очень большой разницей по сравнению с верхним квартилем. Изобразим точки, плотно окруженные объектами OSM, с помощью __folium__."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quant_99 = dist_features.loc[idx_full, 'tagged_points_in_0.25'].quantile(0.99)\n",
    "dense_mask = dist_features.loc[idx_full, 'tagged_points_in_0.25'] > quant_99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dense_zones = zones_df.loc[idx_full][dense_mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#начальным положением зададим Красную площадь\n",
    "fmap = folium.Map([55.753722, 37.620657])\n",
    "\n",
    "for _, row in dense_zones.iterrows():\n",
    "    folium.CircleMarker([row.lat_c, row.lon_c], radius=4).add_to(fmap)\n",
    "    \n",
    "fmap    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Основная масса точек лежит близко к центру Москвы, что, как мне кажется оправдывает такую плотность объектов в участке.\n",
    "\n",
    "__Посмотрим на другие признаки из dist_features__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_features.loc[idx_full].describe().iloc[:, 17:33]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_features.loc[idx_full].describe().iloc[:, 33:49]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_features.loc[idx_full].describe().iloc[:, 49:65]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_features.loc[idx_full].describe().iloc[:, 65:81]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_features.loc[idx_full].describe().iloc[:, 81:97]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_features.loc[idx_full].describe().iloc[:, 97:113]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В целом опять внимание привлекает разница между верхним квартилем и максимальным значением признаков. __Еще раз посмотрим на карту__, на этот раз отметим центры участков, которые плотно окружены магазинами."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quant_99 = dist_features.loc[idx_full, 'shop_points_in_0.25'].quantile(0.99)\n",
    "dense_mask = dist_features.loc[idx_full, 'shop_points_in_0.25'] > quant_99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dense_zones = zones_df.loc[idx_full][dense_mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#начальным положением зададим Красную площадь\n",
    "fmap = folium.Map([55.753722, 37.620657])\n",
    "\n",
    "for _, row in dense_zones.iterrows():\n",
    "    folium.CircleMarker([row.lat_c, row.lon_c], radius=4, color='#DC143C').add_to(fmap)\n",
    "    \n",
    "fmap    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Если вручную потыкать по карте, то можно увидеть, что какие-то отмеченные точки просто находятся рядом с торговыми центрами, отсюда и такое число магазинов и мест с тегом \"food\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__По признакам, связанным с жилыми домами__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_house_features.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Посмотрим на жилые дома на карте__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quant_99 = full_house_features.loc[idx_full, 'n_resident'].quantile(0.99)\n",
    "dense_mask = full_house_features.loc[idx_full, 'n_resident'] > quant_99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dense_zones = zones_df.loc[idx_full][dense_mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#начальным положением зададим Красную площадь\n",
    "fmap = folium.Map([55.753722, 37.620657])\n",
    "\n",
    "for _, row in dense_zones.iterrows():\n",
    "    folium.CircleMarker([row.lat_c, row.lon_c], radius=4, color='#DC143C').add_to(fmap)\n",
    "    \n",
    "fmap    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Большое количество жилых домов особенно не смущает: на карте можно видеть места со спальными районами вроде Балашихи, Химок, Мытищ. Больше 70% данных — это участки без жилых домов вообще, и вот это настораживает больше. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__По числу переходов и БС__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Переходы:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crossing_features.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Здесь аналогичный разрыв между максимумом и верхним квартилем.\n",
    "\n",
    "БС:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_base_stations.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quant_99 = n_base_stations[idx_full].quantile(0.99)\n",
    "dense_mask = n_base_stations[idx_full] > quant_99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dense_zones = zones_df.loc[idx_full][dense_mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#начальным положением зададим Красную площадь\n",
    "fmap = folium.Map([55.753722, 37.620657])\n",
    "\n",
    "for _, row in dense_zones.iterrows():\n",
    "    folium.CircleMarker([row.lat_c, row.lon_c], radius=4, color='#DC143C').add_to(fmap)\n",
    "    \n",
    "fmap    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Снова скопление точек в центре Москвы."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Посмотрим на целевой признак__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Из-за особенностей метрики качества (важен порядок, а не точные значения), можно, вообще говоря, делать одинаковые предсказания на каждый день недели (так делали в Starter Kit) основываясь на среднем числе звонков в день."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calls_train.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Посмотрим на распределение звонков__ по дням недели, по будням/выходным и просто в среднем по дням."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__По дням недели__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax_list = plt.subplots(3, 3, figsize=(15, 10))\n",
    "\n",
    "n = 0\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        sns.kdeplot(calls_train['calls_wd{}'.format(n)], ax=ax_list[i][j]);\n",
    "        n += 1\n",
    "        if n > 6:\n",
    "            break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Уже по графику видно, что распределение ассимметрично. Можно формально проверить на нормальность и ожидаемо отвергнуть гипотезу о нормальности. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import normaltest\n",
    "\n",
    "alpha = 0.05\n",
    "_, pvalue = normaltest(calls_train['calls_wd0'])\n",
    "pvalue > alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__По будням/выходным__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.kdeplot(calls_train['calls_workday'])\n",
    "sns.kdeplot(calls_train['calls_weekend'], color='red');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__В среднем по всем дням__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.kdeplot(calls_train['calls_daily']);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Главный вывод из графиков распределения: целевые признаки не распределены нормально, значения сильно группируются вблизи нуля. Возможно, стоит преобразовать целевой признак, чтобы избавиться от сильной ассимметричности."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Посмотрим на boxplot звонков по дням недели__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.boxplot(data=[calls_train['calls_wd{}'.format(i)] for i in range(7)]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Здесь хорошо видны выбросы: по всем дням недели есть зоны, у которых число звонков заходит сильно выше верхнего квартиля. \n",
    "\n",
    "__Посмотрим поближе__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quant_99 = [calls_train['calls_wd{}'.format(i)].quantile(0.99) for i in range(7)]\n",
    "mx = [calls_train['calls_wd{}'.format(i)].max() for i in range(7)]\n",
    "plt.scatter(x=range(7), y=quant_99, c='red')\n",
    "plt.scatter(x=range(7), y=mx, c='black');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__И на boxplot звонков по будням/выходным__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.boxplot(data=[calls_train.calls_workday, calls_train.calls_weekend]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Из этих графиков можно сделать вывод о том, что есть участки, у которых среднее число вызовов нехарактерно большое. Причем, скорее всего, дело не в самих участках, а в каком-то ЧП, которое там произошло. Пожар, например, из-за которого могло резко увеличиться число звонков. Думаю, что при использовании линейной регрессии такие участки нужно исключить из выборки из-за нетипичности.\n",
    "\n",
    "Еще можно посмотреть на эти участки на карте. Они не локализуются, то есть если это все-таки ЧП, то разные ЧП. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quant_99 = calls_train['calls_daily'].quantile(0.99)\n",
    "mask = calls_train['calls_daily'] > quant_99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "many_calls_zones = zones_df.loc[idx_train][mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#начальным положением зададим Красную площадь\n",
    "fmap = folium.Map([55.753722, 37.620657])\n",
    "\n",
    "for _, row in many_calls_zones.iterrows():\n",
    "    folium.CircleMarker([row.lat_c, row.lon_c], radius=4, color='#DC143C').add_to(fmap)\n",
    "    \n",
    "fmap    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Попробуем посмореть, как коррелируют между собой признаки и как они влияют на целевой__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Для dist_features__\n",
    "\n",
    "tagged_points будут коррелировать со всеми, так что посмотрим на корреляцию признаков вида \"число объектов в заданном радиусе\" без tagged_points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_list = ['railway', 'public_transport', 'entertainment', 'food', 'leisure', 'shop']\n",
    "\n",
    "cmap = sns.cm.rocket_r\n",
    "sns.heatmap(dist_features.loc[idx_full, ['{}_points_in_0.1'.format(p) for p in p_list]]\n",
    "            .corr('kendall'), cmap=cmap);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.heatmap(dist_features.loc[idx_full, ['{}_points_in_0.25'.format(p) for p in p_list]]\n",
    "            .corr(), cmap=cmap);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Тут никаких полезных инсайтов, только бытовые наблюдения: коррелирует транспорт (в городе вокзалы часто рядом с метро и автобусными остановками), самые сильные корреляции — c магазинами (если в участке много магазинов, то там, скорее всего, есть большой торговый центр, в котором есть фуд-корты и которые стоят недалеко от остановок/метро). Возьмем побольше признаков."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Для dist_features, house_features и base_stations__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.heatmap(full_features.loc[idx_full, ['{}_points_in_0.25'.format(p) for p in p_list] +\n",
    "                                        ['n_resident', 'n_old', 'n_ind_proj', 'n_auto_fire', 'n_without_gas'] +\n",
    "                                        ['n_base_stations']]\n",
    "            .corr(), cmap=cmap);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Число жилых домов сильно коррелирует с числом домов старше 1970 года и домов без газоснабжения. Из интересного вижу разве что корреляцию между числом базовых станций и магазинами и местами, где можно поесть: это, возможно, опять можно объяснить торговыми центрами."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Посмотрим на влияние признаков на целевой__\n",
    "\n",
    "Поскольку распределение звонков по дням недели слабо отличается, посмотрим на звонки не по дням недели, а по будням/выходным и в среднем по всем дням."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.heatmap(pd.concat([calls_train[['calls_daily', 'calls_workday', 'calls_weekend']],\n",
    "                       train_features.loc[:, ['{}_points_in_0.25'.format(p) for p in p_list] +\n",
    "                                             ['{}_mean_dist_k_5'.format(p) for p in p_list]] \n",
    "                      ],\n",
    "                       axis=1)\n",
    "                       .corr(), cmap=cmap);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Здесь звонки сильнее всего коррелируют с объектами, относящимися к общественному транспорту, с магазинами и кафе/ресторанами, причем положительно, хотя было предположение, что в случае с магазинами должно быть наоборот. Возможное объяснение: в торговых центрах могут быть пожары, при которых звонят сразу много людей (об этом выше), а рядом с остановками — аварии.\n",
    "\n",
    "При этом со средним расстоянием до разных объектов звонки коррелируют слабо, то есть либо связь сильно нелинейная, либо связь в принципе слабая. Ниже посмотрим на графики; возможно, эти признаки можно будет вообще убрать."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.heatmap(pd.concat([calls_train[['calls_daily', 'calls_workday', 'calls_weekend']],\n",
    "                       train_features.loc[:, ['{}_min'.format(p) for p in p_list] +\n",
    "                                             ['n_crossing'] +\n",
    "                                             ['n_crossing_uncontr'] +\n",
    "                                             ['kremlin_dist']]\n",
    "                      ],\n",
    "                       axis=1)\n",
    "                       .corr(), cmap=cmap);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сильнее всего звонки коррелируют с числом переходов. С минимальным расстоянием до объектов, включая Кремль, по силе корреляция меньше."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.heatmap(pd.concat([calls_train[['calls_daily', 'calls_workday', 'calls_weekend']],\n",
    "                       train_features.loc[:, ['n_resident', 'n_old', 'n_ind_proj',\n",
    "                                              'n_auto_fire', 'n_without_gas', 'n_base_stations']]\n",
    "                      ],\n",
    "                       axis=1)\n",
    "                       .corr(), cmap=cmap);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сильнее всего звонки коррелируют с числом БС. А вот с домами довольно слабо."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.heatmap(pd.concat([calls_train[['calls_daily', 'calls_workday', 'calls_weekend']],\n",
    "                       train_features.loc[:, ['has_clinic', 'has_police', 'has_school']]\n",
    "                      ],\n",
    "                      axis=1)\n",
    "                      .corr(method='spearman'), cmap=cmap);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Есть слабая корреляция звонков с числом больниц."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Посмотрим на графики зависимости__ среднего числа звонков в день от признаков, с которыми имеется наиболее сильная корреляция."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Влияние числа объектов в радиусе 250 м__:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax_list = plt.subplots(2, 3, figsize=(16, 10))\n",
    "\n",
    "ax_num = 0\n",
    "for p in p_list:\n",
    "    sns.regplot(data=pd.concat([train_features, calls_train], axis=1),\n",
    "                x='{}_points_in_0.25'.format(p), \n",
    "                y='calls_daily',\n",
    "                ax=ax_list.reshape(-1, 1)[ax_num][0]);\n",
    "    ax_num += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Что-то похожее на тренд, если не учитывать явные выбросы, видится разве что для shop_points и food_points. У остальных признаков виден сильный разброс относительно регрессионной прямой, которую строит Seaborn."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Влияние числа переходов, БС и жилых домов__:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax_list = plt.subplots(2, 2, figsize=(16, 8))\n",
    "\n",
    "sns.regplot(data=pd.concat([train_features, calls_train], axis=1),\n",
    "            x='n_crossing', \n",
    "            y='calls_daily',\n",
    "            ax=ax_list[0][0]);\n",
    "sns.regplot(data=pd.concat([train_features, calls_train], axis=1),\n",
    "            x='n_crossing_uncontr', \n",
    "            y='calls_daily',\n",
    "            ax=ax_list[0][1]);\n",
    "sns.regplot(data=pd.concat([train_features, calls_train], axis=1),\n",
    "            x='n_base_stations', \n",
    "            y='calls_daily',\n",
    "            ax=ax_list[1][0]);\n",
    "sns.regplot(data=pd.concat([train_features, calls_train], axis=1),\n",
    "            x='n_resident', \n",
    "            y='calls_daily',\n",
    "            ax=ax_list[1][1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "У n_crossing и n_crossing_uncontr разброс относительно прямой и не видно особого тренда, как и у n_resident. Более плотно точки находятся около прямой для n_base_stations, если не учитывать выбросы.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Посмотрим на несколько графиков зависимости__ среднего числа звонков от признаков, с которыми слабая корреляция."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Влияние минимального расстояния до объектов:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax_list = plt.subplots(2, 3, figsize=(16, 10))\n",
    "\n",
    "ax_num = 0\n",
    "for p in p_list:\n",
    "    sns.regplot(data=pd.concat([train_features, calls_train], axis=1),\n",
    "                x='{}_min'.format(p), \n",
    "                y='calls_daily',\n",
    "                ax=ax_list.reshape(-1, 1)[ax_num][0]);\n",
    "    ax_num += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Признаки выглядят неинформативными, связи на графиках не видно. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Влияние среднего расстояния до объектов:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax_list = plt.subplots(2, 3, figsize=(16, 10))\n",
    "\n",
    "ax_num = 0\n",
    "for p in p_list:\n",
    "    sns.regplot(data=pd.concat([train_features, calls_train], axis=1),\n",
    "                x='{}_mean_dist_k_5'.format(p), \n",
    "                y='calls_daily',\n",
    "                ax=ax_list.reshape(-1, 1)[ax_num][0]);\n",
    "    ax_num += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Влияние числа объектов из full_house_info:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax_list = plt.subplots(2, 2, figsize=(16, 8))\n",
    "\n",
    "sns.regplot(data=pd.concat([train_features, calls_train], axis=1),\n",
    "            x='n_old', \n",
    "            y='calls_daily',\n",
    "            ax=ax_list[0][0]);\n",
    "sns.regplot(data=pd.concat([train_features, calls_train], axis=1),\n",
    "            x='n_ind_proj', \n",
    "            y='calls_daily',\n",
    "            ax=ax_list[0][1]);\n",
    "sns.regplot(data=pd.concat([train_features, calls_train], axis=1),\n",
    "            x='n_auto_fire', \n",
    "            y='calls_daily',\n",
    "            ax=ax_list[1][0]);\n",
    "sns.regplot(data=pd.concat([train_features, calls_train], axis=1),\n",
    "            x='n_without_gas', \n",
    "            y='calls_daily',\n",
    "            ax=ax_list[1][1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Стоит как минимум оставить n_ind_proj. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Еще одно наблюдение__\n",
    "\n",
    "Посмотрим на то, как выглядит график зависимости звонков от одного из \"хороших\" признаков, например, от числа магазинов в радиусе 250 м, если из звонков выбросить аномальные примеры:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.regplot(data=pd.concat([train_features,\n",
    "                            calls_train[calls_train.calls_daily < calls_train.calls_daily.quantile(0.995)]], axis=1),\n",
    "                            x='shop_points_in_0.25', \n",
    "                            y='calls_daily');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "И вот здесь становится грустно, потому что основная масса неаномальных точек просто кластеризуется в левом нижнем углу, а оставшиеся точки разбросаны относительно построенной прямой и никакого линейного тренда, да и в принципе тренда, не выделяется. \n",
    "\n",
    "Еще пример. Здесь вроде как можно увидеть тренд (по крайней мере, с ростом числа БС растет минимальное число звонков), но все еще с разбросом:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.regplot(data=pd.concat([train_features,\n",
    "                            calls_train[calls_train.calls_daily < calls_train.calls_daily.quantile(0.995)]], axis=1),\n",
    "                            x='n_base_stations', \n",
    "                            y='calls_daily');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Выводы по анализу:__\n",
    "\n",
    "1) Видны выбросы в числе звонков, что может объясняться крупными, но редкими ЧП. Если использовать линейную модель, то их, думаю, не стоит брать.  \n",
    "\n",
    "2) Про возможные причины корреляций и зависимостей — в промежуточных выводах выше. Субъективно выделяются менее информативные признаки, которые можно проверить на lasso-регрессии, если все-таки использовать в качестве модели регрессию, и отбросить, если их веса занулятся. Либо вообще не трогать при использовании случайного леса, например: возможно, из них можно вытащить что-то полезное, если по-честному смотреть на критерий информативности, а не просто глазами на матрицу корреляций и графики. \n",
    "\n",
    "3) Среди относительно информативных признаков нет чудо-признаков, которые бы хорошо в одиночку объясняли целевую переменную. В лучшем случае можно выделить что-то похожее на зависимось и проверить, как самые информативные признаки будут работать вместе. Это влияет на дальнейший ход задачи: если бы было больше времени, наверное, стоило бы потратить часть его на более серьезный data mining, чем у меня: например, погулять по соцсетям и поанализировать посты с геотегами. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Создание признаков II <a name=\"feat_2\"></a>\n",
    "\n",
    "[К началу](#cont)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "При анализе целевого признака были видны выбросы: области с аномально большим числом звонков. __Введем бинарный признак__  \"есть ли по соседству участок с большим числом звонков\". \"По соседству\" пусть означает, что расстояние между центрами участков не больше 550 м (если рассматривать участки как идеальные квадраты ровно 500х500 м, то расстояние было бы 500 м, еще 50 м накинем сверху из-за возможной неточности). \n",
    "\n",
    "Мотивация: с одной стороны, стоит выкинуть аномальные участки, чтобы не портить модель, чувствительную к выбросам. С другой, хочется как-то сохранить часть информации с этих участков. И даже если нетипичные участки не выкидывать, этот признак, возможно, поможет отделить часть участков не с аномально высоким, но с повышенным числом звонков, основываясь на соседстве с этими нетипичными участками. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def many_calls_neighbour(zone_centre, many_calls_zones=many_calls_zones):\n",
    "    '''\n",
    "    in: array zone_centre (координаты центра участка lat_c, lon_c)\n",
    "    out: 1, если участок рядом с участком с высоким числом звонков, 0 иначе\n",
    "    '''\n",
    "    for many_calls_zone in many_calls_zones[['lat_c', 'lon_c']].values:\n",
    "        if distance(zone_centre, many_calls_zone) <= 0.55:\n",
    "            return 1\n",
    "    return 0    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mcn = zones_df[(zones_df.is_target == 1) | (zones_df.is_test == 1)]. \\\n",
    "      apply(lambda row: many_calls_neighbour(np.array([row.lat_c, row.lon_c])), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_features['mcn'] = mcn[idx_train]\n",
    "test_features['mcn'] = mcn[idx_test]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Выбор моделей и метрики <a name=\"model\"></a>\n",
    "\n",
    "[К началу](#cont)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для оценки предсказаний используется тау Кендалла, и в задаче правильное упорядочивание важнее, чем точные численные ответы. Есть предположение о том, что похожие участки со схожими условиями проживания будут иметь похожее число звонков, а непохожесть наоборот позволит сравнивать число звонков в терминах меньше/больше — поэтому попробуем использовать метод ближайших соседей. Еще попробуем применить случайный лес из тех же соображений: лес сможет разбить объекты выборки на группы со схожими параметрами и схожим числом звонков, при этом обе эти модели смогут смоделировать нелинейную связь между признаками и числом звонков. А в конце посмотрим на то, как заработает градиентный бустинг. \n",
    "\n",
    "Метрика качества будет, понятно, тау Кендалла, раз она используется при оценке предсказаний.\n",
    "\n",
    "Так как предсказания по дням недели оцениваются независимо, будем пока в качестве целевого признака брать среднее число звонков по всем дням (calls_daily). То есть получится, что для каждого участка предсказанное число звонков с понедельника по воскресенье будет одинаковым. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Предобработка данных и обучение <a name=\"learning\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### KNeighborsRegressor <a name=\"knr\"></a>\n",
    "\n",
    "[К началу](#cont)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "from sklearn.model_selection import validation_curve, learning_curve, KFold\n",
    "from sklearn.pipeline import Pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В признаках вида prt_* много NaN (больше половины), заменять их средним/медианным значением или 0 будет слишком неинформативно, поэтому их в качестве признаков брать не будем. Бинарные признаки, слабо коррелирующие с целевой переменной, тоже попробуем не использовать.\n",
    "\n",
    "Метод ближайших соседей чувствителен к масштабированию, так что используем StandardScaler. Участки с аномально большим числом звонков тоже пока выбросим. От большой размерности признакового пространства KNN может пострадать, так что используем PCA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Построим кривую валидации и подберем параметр n_neighbors__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kendalltau_scorer(estimator, X, y):\n",
    "    '''\n",
    "    Возвращает коэффициент корреляции Кендалла;\n",
    "    такой формат используется для совместимости с validation_curve и learning_curve\n",
    "    '''\n",
    "    from scipy.stats import kendalltau\n",
    "    return kendalltau(estimator.predict(X), y).correlation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Подбирать будем только один гиперпараметр модели: число ближайших соседей, по которым путем усреднения числа звонков число звонков будет предсказываться для нового объекта. Подбирать будем на кросс-валидации (пятикратной, потому что данных не слишком много; с перемешиванием, раз временных меток у данных нет и сохранять временную структуру не надо; с помощью простого перебора). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "drop_cols = ['prt_old', 'prt_ind_proj', 'prt_without_gas', 'prt_auto_fire',\n",
    "             'prt_crossing', 'has_clinic', 'has_school', 'has_police']\n",
    "knr_params = range(1, 61)\n",
    "knr_pipe = Pipeline([('scaler', StandardScaler()), ('pca', PCA(n_components=0.95, random_state=27)),\n",
    "                     ('knr', KNeighborsRegressor())])\n",
    "\n",
    "train_score, test_score = validation_curve(knr_pipe, \n",
    "                                           train_features.drop(drop_cols, axis=1).drop(many_calls_zones.index),\n",
    "                                           calls_train['calls_daily'].drop(many_calls_zones.index),\n",
    "                                           'knr__n_neighbors', knr_params,\n",
    "                                           cv=KFold(n_splits=5, shuffle=True, random_state=27), \n",
    "                                           scoring=kendalltau_scorer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(knr_params, train_score.mean(axis=1), label='train score')\n",
    "plt.plot(knr_params, test_score.mean(axis=1), label='test score');\n",
    "plt.xlabel('n_neighbors')\n",
    "plt.xticks(np.arange(2, 61, 6))\n",
    "plt.ylabel('kendall tau')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В районе n_neighbors=8 у среднего скора на тесте максимум, кривая для теста близко к кривой для трейна, значит, можно предположить, что модель не переобучилась. Возьмем n_neighbors, отвечающий за максимальный средний тестовый скор."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_neighbors = knr_params[np.argmax(test_score.mean(axis=1))]\n",
    "n_neighbors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "И посмотрим на максимальный средний тестовый скор."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_score.mean(axis=1).max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Построим кривую обучения с выбранным n_neighbors__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "train_sizes=np.linspace(0.02, 1, 20)\n",
    "knr_pipe = Pipeline([('scaler', StandardScaler()), ('pca', PCA(n_components=0.95, random_state=27)),\n",
    "                     ('knr', KNeighborsRegressor(n_neighbors=n_neighbors))])\n",
    "N_train, train_score, test_score = learning_curve(knr_pipe, \n",
    "                                         train_features.drop(drop_cols, axis=1).drop(many_calls_zones.index),\n",
    "                                         calls_train['calls_daily'].drop(many_calls_zones.index),\n",
    "                                         train_sizes=train_sizes,\n",
    "                                         cv=KFold(n_splits=5, shuffle=True, random_state=27), \n",
    "                                         scoring=kendalltau_scorer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(N_train, train_score.mean(axis=1), label='train score')\n",
    "plt.plot(N_train, test_score.mean(axis=1), label='test score');\n",
    "plt.xlabel('train size')\n",
    "plt.ylabel('kendall tau')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Кривые для трейна и теста еще не сошлись, у обоих видна тенденция на рост с увеличением размера выборки. \n",
    "\n",
    "Попробуем __сделать предсказания этой модели на тесте и обучить еще одну__ на всех данных из обучающей выборки (до этого использовались только целевые участки). Признаков, соответственно, станет меньше, потому что некоторые получать долго и дорого и они считались только для целевых участков.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_submission_file(y_pred, filename, idx_test=idx_test):\n",
    "    \n",
    "    '''\n",
    "    Записывает предсказанные результаты в файл;\n",
    "    y_pred -- предсказания для числа звонков в среднем по всем дням\n",
    "    '''\n",
    "    from collections import OrderedDict\n",
    "    \n",
    "    calls = ['calls_wd{0}'.format(i) for i in range(0, 7)]\n",
    "    pred_df = pd.DataFrame(OrderedDict([(col_name, y_pred) for col_name in calls]), index=idx_test)\n",
    "    pred_df.to_csv(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Масштабирование, PCA__\n",
    "\n",
    "В пайплайне выше масштабировался и бинарный признак mcn, потому что не хотелось с этим возиться, здесь же его масштабировать не будем. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = StandardScaler()\n",
    "\n",
    "drop_cols = ['prt_old', 'prt_ind_proj', 'prt_without_gas', 'prt_auto_fire',\n",
    "             'prt_crossing', 'has_clinic', 'has_school', 'has_police', 'mcn']\n",
    "\n",
    "scaler.fit(train_features.drop(drop_cols, axis=1).drop(many_calls_zones.index))\n",
    "train_features_scaled = np.concatenate([scaler.transform(train_features.drop(drop_cols, axis=1).drop(many_calls_zones.index)), \n",
    "                                       train_features['mcn'].drop(many_calls_zones.index).values.reshape(-1, 1)], axis=1)\n",
    "test_features_scaled = np.concatenate([scaler.transform(test_features.drop(drop_cols, axis=1)), \n",
    "                                       test_features['mcn'].values.reshape(-1, 1)], axis=1)    \n",
    "\n",
    "pca = PCA(n_components=0.95, random_state=27)\n",
    "\n",
    "pca.fit(train_features_scaled)\n",
    "train_features_pca = pca.transform(train_features_scaled)\n",
    "test_features_pca = pca.transform(test_features_scaled)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Обучение с лучшим гиперпараметром, предсказание для теста__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "knr = KNeighborsRegressor(n_neighbors=n_neighbors)\n",
    "\n",
    "knr.fit(train_features_pca, calls_train['calls_daily'].drop(many_calls_zones.index))\n",
    "calls_pred = knr.predict(test_features_pca)\n",
    "write_submission_file(calls_pred, 'knr.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Результат на LB__\n",
    "\n",
    "34 место из 63 на тот момент.\n",
    "\n",
    "<img src=\"https://habrastorage.org/webt/tx/ye/1o/txye1o9annxuglvpxoh-goa6u_i.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### KNeighborsRegressor на расширенной выборке <a name=\"knr_ext\"></a>\n",
    "\n",
    "[К началу](#cont)\n",
    "\n",
    "__Дополним два признака и создадим новые выборки__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "n_base_stations_non_target = zones_df[zones_df.is_target == 0] \\\n",
    "                            .apply(lambda row: obj_count((row.lat_bl, row.lon_bl, row.lat_tr, row.lon_tr),\n",
    "                                                stat_coords), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_base_stations_ext = pd.concat([n_base_stations, n_base_stations_non_target]).sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_train_ext = zones_df[zones_df.is_test == 0].index\n",
    "calls_train_ext = zones_df.loc[idx_train_ext, :].iloc[:, 8:]\n",
    "\n",
    "quant_99 = calls_train_ext['calls_daily'].quantile(0.99)\n",
    "mask = calls_train_ext['calls_daily'] > quant_99\n",
    "\n",
    "many_calls_zones_ext = zones_df.loc[idx_train][mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "mcn_non_target = zones_df[zones_df.is_target == 0]. \\\n",
    "                 apply(lambda row: many_calls_neighbour(np.array([row.lat_c, row.lon_c]),\n",
    "                                                        many_calls_zones=many_calls_zones_ext), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mcn_ext = pd.concat([mcn, mcn_non_target]).sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_features_ext = pd.concat([dist_features,\n",
    "                               crossing_features,\n",
    "                               n_base_stations_ext.rename('n_base_stations'),\n",
    "                               mcn_ext.rename('mcn')],\n",
    "                               axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_features_ext = full_features_ext.loc[idx_train_ext]\n",
    "test_features_ext = full_features_ext.loc[idx_test]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Посмотрим на кривую валидации__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "knr_params = range(1, 21)\n",
    "knr_pipe = Pipeline([('scaler', StandardScaler()), ('pca', PCA(n_components=0.95, random_state=27)),\n",
    "                     ('knr', KNeighborsRegressor())])\n",
    "\n",
    "train_score, test_score = validation_curve(knr_pipe, \n",
    "                                           train_features_ext.drop('prt_crossing', axis=1).drop(many_calls_zones_ext.index),\n",
    "                                           calls_train_ext['calls_daily'].drop(many_calls_zones_ext.index),\n",
    "                                           'knr__n_neighbors', knr_params,\n",
    "                                           cv=KFold(n_splits=5, shuffle=True, random_state=27), \n",
    "                                           scoring=kendalltau_scorer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(knr_params, train_score.mean(axis=1), label='train score')\n",
    "plt.plot(knr_params, test_score.mean(axis=1), label='test score');\n",
    "plt.xlabel('n_neighbors')\n",
    "plt.xticks(np.arange(2, 21))\n",
    "plt.ylabel('kendall tau')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Кривые не сходятся, как в прошлый раз, но не сильно отличаются."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_neighbors = knr_params[np.argmax(test_score.mean(axis=1))]\n",
    "n_neighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_score.mean(axis=1).max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Масштабирование, PCA__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = StandardScaler()\n",
    "\n",
    "scaler.fit(train_features_ext.drop('prt_crossing', axis=1).drop(many_calls_zones_ext.index))\n",
    "train_features_ext_scaled = np.concatenate([scaler.transform(train_features_ext.drop('prt_crossing', axis=1).drop(many_calls_zones_ext.index)), \n",
    "                            train_features_ext['mcn'].drop(many_calls_zones_ext.index).values.reshape(-1, 1)], axis=1)\n",
    "test_features_ext_scaled = np.concatenate([scaler.transform(test_features_ext.drop('prt_crossing', axis=1)), \n",
    "                           test_features_ext['mcn'].values.reshape(-1, 1)], axis=1)    \n",
    "\n",
    "pca = PCA(n_components=0.95, random_state=27)\n",
    "\n",
    "pca.fit(train_features_ext_scaled)\n",
    "train_features_ext_pca = pca.transform(train_features_ext_scaled)\n",
    "test_features_ext_pca = pca.transform(test_features_ext_scaled)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Обучение с лучшим гиперпараметром, предсказание для теста__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "knr = KNeighborsRegressor(n_neighbors=n_neighbors)\n",
    "\n",
    "knr.fit(train_features_ext_pca, calls_train_ext['calls_daily'].drop(many_calls_zones_ext.index))\n",
    "calls_pred = knr.predict(test_features_ext_pca)\n",
    "write_submission_file(calls_pred, 'knr_extended.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Результаты на LB__\n",
    "\n",
    "Все здорово просело. По всей видимости, учитывать нецелевые участки все-таки не стоит, либо стоит как-то по-хитрому, например, отсекая участки с совсем маленьким числом звонков и брать топ-n по звонкам. Если останется время, то можно попробовать и это. \n",
    "\n",
    "<img src=\"https://habrastorage.org/webt/id/6e/-u/id6e-urdmvodkr7_db8vxsfdjgq.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RandomForestRegressor <a name=\"rfr\"></a>\n",
    "\n",
    "[К началу](#cont)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Здесь строить кривые валидации и обучения уже не будем, просто посмотрим на cross_val_score, параметры у валидации те же. Гиперпараметры для леса тоже особо подбирать не будем и воспользуемся советами из [статьи курса](https://habrahabr.ru/company/ods/blog/324402/#2-sluchaynyy-les) насчет размера подмножества признаков и минимального числа объектов в листах. Деревьев попробуем взять 500 (лес не переобучается от количества деревьев, так что можно и больше).\n",
    "\n",
    "От выбросов лес не пострадает, так что и участки с аномальным числом звонков выбрасывать не будем."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drop_cols = ['prt_old', 'prt_ind_proj', 'prt_without_gas', 'prt_auto_fire', 'prt_crossing']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "n_features = train_features.drop(drop_cols, axis=1).shape[1]\n",
    "rfr = RandomForestRegressor(n_estimators=500, min_samples_leaf=5, max_features=round(n_features / 3), random_state=27)\n",
    "cvs = cross_val_score(estimator=rfr, X=train_features.drop(drop_cols, axis=1), y=calls_train['calls_daily'],\n",
    "                      scoring=kendalltau_scorer, cv=KFold(n_splits=5, shuffle=True, random_state=27))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvs.mean(), cvs.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rfr.fit(train_features.drop(drop_cols, axis=1), calls_train['calls_daily'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Посмотрим на значимость признаков__\n",
    "\n",
    "Некоторые новые признаки, которые создавались вручную не на основе StarterKit, оказались в топ-20: число зон-соседей с большим числом звонков, число жилых домов в участке, число БС в участке. Еще для леса оказались важными признаки вроде минимальное/максимальное расстояние до какого-то объекта, хотя предполагалось, что эти признаки заработают хуже, чем число каких-то объектов в участке. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "importances = rfr.feature_importances_\n",
    "indices = np.argsort(importances)[::-1]\n",
    "\n",
    "sns.barplot(y=train_features.drop(drop_cols, axis=1).columns[indices[:21]],\n",
    "            x=importances[indices[:21]],\n",
    "            orient='h');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calls_pred = rfr.predict(test_features.drop(drop_cols, axis=1))\n",
    "write_submission_file(calls_pred, 'rfr.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Результаты на LB__\n",
    "\n",
    "14 место из 65 на тот момент.\n",
    "\n",
    "<img src=\"https://habrastorage.org/webt/jk/-h/8u/jk-h8uembrj7xszlnhswynsjdbk.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### XGBRegressor <a name=\"xgb\"></a>\n",
    "\n",
    "[К началу](#cont)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install xgboost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xgboost.sklearn import XGBRegressor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В качестве базовых алгоритмов возьмем деревья. Настраивать гиперпараметры будем в два этапа: сначала для базовых алгоритмов (max_depth, максимальную глубину деревьев; colsample_bytree, долю признаков, на которой будет обучаться каждое дерево, и gamma,  минимальное уменьшение функции потерь, при котором вершина дерева еще будет разбиваться). \n",
    "\n",
    "max_depth и gamma напрямую контролируют сложность базовой модели, а colsample_bytree добавляет шума для большей устойчивости модели ([документация XGBoost](http://xgboost.readthedocs.io/en/latest/how_to/param_tuning.html?highlight=parameters)). Интервалы перебора max_depth и colsample_bytree берутся навскидку, gamma с предварительной оценкой.\n",
    "\n",
    "Потом настроим темп обучения и число деревьев в модели. Также оставим отложенную выборку и результаты обоих этапов настройки гиперпараметров оценим на ней."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Оценим возможный интервал для gamma__\n",
    "\n",
    "__Подберем первый блок параметров__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Решается задача регрессии, в качестве функции потерь будет по умолчанию использоваться RMSE. \n",
    "\n",
    "$$RMSE = \\sqrt{\\frac{1}{n}\\Sigma_{i=1}^{n}{({\\hat{y_i} - y_i})^2}}$$ \n",
    "\n",
    "Поскольку предсказываются доли, то максимальной RMSE будет, когда всегда предсказывается число звонков, равное единице, когда на самом деле звонков всегда ноль. Тогда RMSE лежит в интервале от 0 до 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "X_train, X_valid, y_train, y_valid = train_test_split(train_features.drop(drop_cols, axis=1),\n",
    "                                                      calls_train['calls_daily'],\n",
    "                                                      test_size=0.2,\n",
    "                                                      random_state=27)\n",
    "\n",
    "xgbr_params = {'max_depth': np.arange(3, 16, 3),\n",
    "               'colsample_bytree': [0.3, 0.6, 0.9], \n",
    "               'gamma': np.append(np.logspace(-4, -1, 4), 0)}\n",
    "\n",
    "cv = GridSearchCV(estimator=XGBRegressor(random_state=27), param_grid=xgbr_params, scoring=kendalltau_scorer,\n",
    "                  cv=KFold(n_splits=5, shuffle=True, random_state=27), n_jobs=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "cv.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv.best_params_, cv.best_score_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Оценим на отложенной выборке__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kendalltau(cv.predict(X_valid), y_valid).correlation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Посмотрим на второй блок параметров__\n",
    "\n",
    "Попробуем увеличить число деревьев и снизить темп обучения (по умолчанию 0.1), learning_rate, при этом увеличив число итераций, num_boost_round (по умолчанию 10). Будем использовать оптимальные параметры с кросс-валидации выше."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgbr_params = {'n_estimators': [200, 300, 400, 500],\n",
    "               'learning_rate': [0.01, 0.001]}\n",
    "set_params = {'max_depth': 15, \n",
    "              'gamma': 0,\n",
    "              'colsample_bytree': 0.6,\n",
    "              'random_state': 27,\n",
    "              'num_boost_round': 20}\n",
    "\n",
    "cv = GridSearchCV(estimator=XGBRegressor(**set_params),\n",
    "                  param_grid=xgbr_params, scoring=kendalltau_scorer,\n",
    "                  cv=KFold(n_splits=5, shuffle=True, random_state=27), n_jobs=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "cv.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv.best_params_, cv.best_score_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__И оценим на отложенной выборке__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kendalltau(cv.predict(X_valid), y_valid).correlation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Обучим xgb на всей выборке и сделаем предсказания на тестовой__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_params = {'max_depth': 15, \n",
    "                'gamma': 0,\n",
    "                'colsample_bytree': 0.6,\n",
    "                'random_state': 27,\n",
    "                'num_boost_round': 20,\n",
    "                'learning_rate': 0.01,\n",
    "                'n_estimators': 500}\n",
    "xgbr = XGBRegressor(**final_params)\n",
    "xgbr.fit(train_features.drop(drop_cols, axis=1), calls_train['calls_daily'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "write_submission_file(xgbr.predict(test_features.drop(drop_cols, axis=1)), 'xgbr.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Результат на LB__\n",
    "\n",
    "Хуже, чем у леса. Вероятно, стоит еще повозиться с гиперпараметрами, например, подбирать число итераций и параметры ранней остановки тоже на кросс-валидации. \n",
    "\n",
    "<img src=\"https://habrastorage.org/webt/kz/cz/af/kzczafhzo7irgs2rmc7_9lt4a_i.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Почти выводы <a name=\"alm_conc\"></a>\n",
    "\n",
    "[К началу](#cont)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Наткнулась на [публичный кернел](http://nbviewer.jupyter.org/github/VlasovKirill/mts-geohack/blob/master/2-public_ranking.ipynb) от Кирилла Власова и воспользовалась его идеей преобразовать целевой признак так, чтобы его распределение было менее ассимметричное и больше похоже на нормальное, что, в предположении, должно лучше сохранять порядок зон по числу звонков.\n",
    "\n",
    "Преобразование — корень от логарифма (log1p) целевого признака. На графики, где видно, как изменяется целевой признак, можно посмотреть в кернеле. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "n_features = train_features.drop(drop_cols, axis=1).shape[1]\n",
    "rfr = RandomForestRegressor(n_estimators=500, min_samples_leaf=5, max_features=round(n_features / 3), random_state=27)\n",
    "cvs = cross_val_score(estimator=rfr, X=train_features.drop(drop_cols, axis=1), y=np.sqrt(np.log1p(calls_train['calls_daily'])),\n",
    "                      scoring=kendalltau_scorer, cv=KFold(n_splits=5, shuffle=True, random_state=27))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvs, cvs.mean(), cvs.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rfr.fit(train_features.drop(drop_cols, axis=1), np.sqrt(np.log1p(calls_train['calls_daily'])))\n",
    "calls_pred = rfr.predict(test_features.drop(drop_cols, axis=1))\n",
    "calls_pred = np.expm1(calls_pred ** 2)\n",
    "write_submission_file(calls_pred, 'rfr_calls_modelled.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Результаты на LB__\n",
    "\n",
    "С лесом и преобразованием получилось 10 место из 71 на тот момент.\n",
    "\n",
    "<img src=\"https://habrastorage.org/webt/fv/pc/as/fvpcasemde11qmwshqq9rddp1ly.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Текущий рейтинг__\n",
    "\n",
    "Посылок больше, чем тех, которые выносились в этот проект: в проекте пропускались первые пробные посылки на основе стартера."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://habrastorage.org/webt/k2/nt/rk/k2ntrklewoe82f9xhc8hlkkzme4.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Выводы <a name=\"conc\"></a>\n",
    "\n",
    "[К началу](#cont)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Если по-честному, то корреляция в 67% выглядит довольно маленькой и ненадежной, если решение всерьез применять на практике. Наверное, отчасти это объясняется сложностью моделируемого процесса, в котором много случайности, хотя есть и какие-то закономерности, которые модель улавливает. Думаю, чтобы решение действительно заработало, нужно побольше данных в обучающей выборке и больше информативных признаков. С признаками вообще получается одновременно и сложнее, и проще, чем в соревнованиях, где внешние данные не разрешены: проще, потому что можно брать вообще все, что одобрят организаторы, и сложнее, потому что нужно много времени тратить на обдумывание того, что нужно искать, сам поиск и предобработку. Зато интересно. \n",
    "\n",
    "Много чего можно улучшить в самом решении, конечно:\n",
    "\n",
    "1) Использовать данные по ЧП (статистику пожаров и ДТП, например). Проблема: сходу такие данные в удобном формате не ищутся, много возни с обработкой. Плюс такие данные полезны для соревнования, но для реального применения на практике — не очень, потому что признаки на их основе были бы по сути построены на статистике из будущего.\n",
    "\n",
    "2) Использовать данные из соцсетей и, возможно, часть нецелевых участков, об этом писала выше.\n",
    "\n",
    "3) Использовать данные, которые хотелось использовать, если бы они нашлись: хотела построить признаки на основе данных по жилым домам с камерами видеонаблюдения во дворе/в подъезде, но нашла такие данные только для Москвы, но не для МО.\n",
    "\n",
    "4) Использовать более сложные модели; получше настроить xgboost, который должен зайти в этой задаче, потому что признаков немного; либо вообще подойти к задаче не как к задаче регрессии ([вот здесь](http://nbviewer.jupyter.org/github/VlasovKirill/mts-geohack/blob/master/2-public_ranking.ipynb) Кирилл Власов подходит к задаче как к задаче ранжирования: соответственно, меняется функция потерь и можно более явно оптимизировать именно порядок, а не точные значения; я о таком подходе не думала.\n",
    "\n",
    "5) Попробовать решать задачу как классификацию: кластеризовать участки из трейна по какому-то подмножеству информативных признаков, таким образом получить кластеры-классы.\n",
    "\n",
    "6) Поработать с целевым признаком: например, отдельно предсказывать разные дни недели, учитывать рабочие/нерабочие дни.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
